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ABSTRACT 

We study numerical convergence in local two-dimensional hydrodynamical simulations 
of self-gravitating accretion discs with a simple cooling law. It is well-known that there 
exists a steady gravito-turbulent state, in which cooling is balanced by dissipation of 
weak shocks, with a net outward transport of angular momentum. Previous results 
indicated that if cooling is too fast (typical time scale 3 , where J7 is the local 
angular velocity), this steady state can not be maintained and the disc will fragment 
into gravitationally bound clumps. We show that, in the two-dimensional local ap- 
proximation, this result is in fact not converged with respect to numerical resolution 
and longer time integration. Irrespective of the cooling time scale, gravito-turbulence 
consists of density waves as well as transient clumps. These clumps will contract be- 
cause of the imposed cooling, and collapse into bound objects if they can survive for 
long enough. Since heating by shocks is very local, the destruction of clumps is a 
stochastic process. High numerical resolution and long integration times are needed 
to capture this behaviour. We have observed fragmentation for cooling times up to 20 
almost a factor 7 higher than in previous simulations. Fully three-dimensional 
simulations with a more realistic cooling prescription are necessary to determine the 
effects of the use of the two-dimensional approximation and a simple cooling law. 

Key words: planets and satellites: formation -planetary systems: protoplanetary 
discs - accretion discs - hydrodynamics - instabilities 



1 INTRODUCTION 

The notion that gaseous planets may form by gravitational 
instability (GI) dates back to iKuiper ( 1951 1 and Cameron 



( 1978[ ). Although revived by|Bossj(jl997 



the disc instability 

model suffered a blow when it was realised that very efficient 
cooling was needed for it to operate (Gammie 20011. This 



makes GI ineffective in planet forming regions (1-20 AU) 
of typical protoplanetary discs (e.g. 'Rafikov '2005 ; Matzner' 
fc Levin]|2005 1. The recent discovery of planets orbiting at 



very large distances from their central star ( Kalas et al. 2008 
[Lafreniere et al.]|2010[ ), where cooling time scales are suffi- 
ciently short, has sparked new interest in the possibility of 
planet formation via disc instability (e.g. Boloy 2009). For 
a review of gravitational instability and its applications in 



planet formation theory, see Durisen et al. ( 2007 1 



A razor-thin self-gravitating disc is unstable to axisym- 
metric perturbations when ( |Safronov||1960| |Toomre||1964[ ): 



Q 



ttGE 



< 1, 



(1) 
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where Cs is the sound speed, k is the epicyclic frequency, 
G is the gravitational constant and E is the surface density. 
For Keplerian discs, k equals the angular velocity of the disc 
Q.. Discs are in fact unstable to non-axisymmetric perturba- 
tions at slightly higher values of Q ( [Papaloizou &: Savonije| 
19911. Dissipation of the resulting spiral waves provides a 



source of heating, which can balance radiative cooling. An 
equilibrium can then be set up ( |Paczynski||1978| , in which 
gravitational and Reynolds stresses are such that the energy 
losses due to cooling are compensated for by dissipation. If 
we parametrise the cooling time scale as 



^cool — , 



(2) 



then in equilibrium we must have that ( |Pringle|1981| |Gam- 
mlel|2001[ ): 



97(7- l)/3' 



(3) 



where a is the usual stress parametrisation ( [Shakura et al.| 
1978 1 and 7 is the ratio of specific heats. This equilibrium of 



gravito-turbulence can therefore lead to significant transport 
of angular momentum in regions of the disc where cooling 
is efficient and Q ~ 1. 



2 S.-J. Paardekooper 



The above parametrisation assumes angular momentum 
transport is local, i.e. that the spiral waves damp close to 
their location of excitation ( Cossins et al.|2d09 Forgan et al 
[2011 j). In principle, there can be a non-local part in the trans- 



port (Balbus & Papaloizou 19991, but Smoothed Particle 



Hydrodynamics (SPH) simulations find that transport is es 
sentially local for low disc masses Mdisc/Af* < 0.25 (Lodato 



& Rice 2004, 2005 ). In the present work, we restrict ourselves 



to local simulations, for which the boundary conditions en- 
sure local transport ( Balbus fc Papaloizou|1999 ). 

Naturally, the critical value for tcooi, which can be trans- 
lated into a critical value of (3, /3c, through equation ([2|, 
below which fragmentation is inevitable, has received much 
attention, since this defines the region in discs where planets 
can form by GI. Gammie ( 2001 1 found that /3c — 3 for 2D 



local simulations with 7 = 2. Rice et al. ( 2005[ ) interpreted 
the critical value for /3 as a maximum value of a the disc can 
sustain (see equation S). They found Omax = 0.06, which 



is in agreement with the /3c found by Gammie (2001 1. Note 



that this implies that /3c depends on 7. 

The simple parametrisation of the cooling time scale 
in terms of the local dynamical time scale is of course an 
oversimplification. There have been several attempts to im- 
prove on the simple cooling law of equation ([2|. [Nelson et al 



( |2000[ ) discuss an implementation of photospheric cooling. 



while the simulations of 'Boss (2001 1 and |Boley et al. ([2006 1 
include full radiative transfer in their global models. |John-| 
[son &: Gammie| ( |2003[ ) pointed out that care has to be taken 
in applying the simple cooling criterion /3 < /3c to discs with 
realistic cooling. The development of non-linear structures 
in the disc can cause the cooling time scale to change signif- 
icantly from that of the initial state. More recently, [Kratter] 
fc Murray-Clay| ( |2011| l and |Rice et"al] ( |2011[ ) studied the ef- 
fect of irradiation by the central star on disc fragmentation. 
It is generally agreed that for these more complicated mod- 
els, similar to the case of a simple cooling law, fast cooling 
is needed to trigger fragmentation. 

Recently, numerical convergence of the critical cooling 
time scale was questioned in Meru & Bate (20111, who 



found fragmentation at much higher values of /3 than pre- 
viously found at lower resolution. It was pointed out in 



Paardekooper et al. (20111 that this behaviour can at least 



partly be explained by the effect of starting from smooth ini- 
tial conditions. In global simulations with constant /3, where 
the cooling time scale therefore depends on radius through 
Q in equation ([2|, the inner parts of the disc get turbulent 
before the outer parts. The edge between the turbulent and 
non-turbulent region of the disc can trigger fragmentation at 
longer cooling times. In local simulations, both /3 and tcooi 
are constant, and therefore no edges appear. However, even 
in this case one has to worry about initial conditions: if the 
time scale to set up gravito-turbulence is longer than the 
cooling time scale, the disc will quickly cool down to Q < 1, 
which triggers the strong linear instability, possibly leading 
again to artificial fragmentation. 

The infiuence of artificial viscosity on fragmentation in 



SPH simulations was discussed in Lodato & Clarke (20111, 



who suggested that artificial heating could play a major role 
in fragmentation even if it constitutes only 5% of the total 



heating. Pickett & Durisen (20071 showed that in locally 



isothermal, grid-based three-dimensional simulations, artifi- 
cial viscosity affects the survival of clumps. It is clear that 



disc fragmentation is a difficult problem to handle numer- 
ically, and therefore it makes sense to go back to the sim- 
plest possible set-up in which fragmentation can be studied, 
namely, the two-dimensional shearing sheet. 

Numerical convergence in local simulations with respect 
to f3c has not been thoroughly discussed so far. [Gammie] 
( |2001[ ) showed convergence with respect to the value of a 
measured in the simulations, but it remains to be seen if the 
fragmentation criterion is independent of numerical resolu- 
tion. This is the subject of the present work. The paper is 
organised as follows: in Sect. [2] we present the basic equa- 
tions in the framework of the shearing sheet, and outline 
the numerical method used to solve these equations in Sect. 
[3] We discuss the performance of the numerical method on 
two test problems in Sect. |4] In Sect. [5] we discuss the ini- 
tial conditions used, and present the results in Sect. [6] We 
discuss the results in Sect. [T] and conclude in Sect.|8l 



2 BASIC EQUATIONS 



Following [Gammie] (|2001[), we use the shearing sheet ap- 
proximation (Goldreich & Lynden-Bell 19651, in which we 



follow a patch of the disc centred at cylindrical coordinates 
(Rjip) = {Rojifio + ^t), where n is the angular velocity at 
R — Ro- On this patch, we define a Cartesian coordinate 
frame defined hy x = R — Ro and y = Ro{ip — ipo — Qt), 
and expand the equations of motion to first order in |x|/_Ro. 
This leads to a linear shear in the patch. The basic equations 
then are the continuity equation: 



9E 
dt 



V ■ Ev = 0, 



(4) 



with E the surface density and v — {vx,Vy)'^ the 2D velocity 
vector, Euler's equation: 

^-Fv-Vv+^ =2gn^sx-2!:!x v-V<E>, (5) 

where p is the 2D pressure, q is the shear parameter {q = 
3/2 in a Keplerian disc), f2 is the angular velocity of the 
coordinate frame, and <I> is the self-gravity potential. Finally, 
we have the equation for internal energy density: 



9e 
dt 



-I- V ■ ev = — pV • V 



en 



(6) 



where the last term represents our simple cooling law. We 
will assume a perfect gas throughout so that p = (7 — l)e. 

The self-gravity potential of a three-dimensional density 
distribution p can be found from 



$(x) 



-G 



p(x')c/^x' 



Taking p — E(x, y)S{z — S), we get 

E(a;', y')dxdy 



(7) 



^(x,y,z) = -G , .g. 

J ^{x - x'Y + {y- y'Y + (2 - 5)2 

from which we see that $(2: = 0) gives the potential due to a 
surface density E, but smoothed over a length 5. Such a po- 
tential would arise when integrating the three-dimensional 
equations in the vertical direction over the thickness of the 
disc. 

For a surface density that is periodic in all directions 
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and for a single Fourier component k : 
equation has the solution 



-27rGf^cxp(-|k|(5). 



{kx, ky) , the above 



(9) 



The smoothing affects Fourier components on scales smaller 
than 5. Note that both |Gamrnie| ( |200ip and |Rice et"aL] ( |2011[ ) 
use 5 = 0, allowing for small scales in the potential (up to 
the grid scale). 

We assume a background state that has constant surface 
density Eo and sound speed Cs,o- This defines both a Q-value 
for the patch (equation ([T])) and an approximate measure of 
the disc thickness H — Cs,o/f2- Note that in the presence of 
self-gravity, H is no longer exactly the pressure scale height. 

It is convenient to split-off the Keplerian part of the 
velocity and work with the perturbed velocity u — {vx , Vy + 



9E ^ as „ 



du 



Su = 0, 



^ 9u „ Vp 
--g^7x- + u.Vu+^ 



(10) 



:gr2u^y-2Sl X u- V*, (11) 



de de 
dt dy 



V • eu = — pV • u 



(12) 



Expressing these in conservation form, defining momenta 
m = TiVx and n — 'S{vy + qVlx), we get: 

9E 



^ 9E ^ dm ^ dn 
dt dy dx dy 



0, 



dm ^ dm d f m 



d_ 

dy 



mn 



dn 



d 



qQx^^ + 
dt dy dx 



de Q '^^ 

dt dy dx 



fmn\ 



d^ 
dy 



2D,n - E 



9$ 

dx ' 



(q~2)nm-T, 



+ P j = 
9$ 



dy 



qfl- 



dx 



(13) 



(14) 



(15) 



(16) 



dy (7-1)/?' 

where e = (m^ -I- n^)/2E -I- e is the total energy density as- 
sociated with the velocity perturbations. Note that we have 
removed the contribution from the shear from the total en- 
ergy. Together with taking no background gradients, this 
ensures that all quantities E, m, n, p, <I> and e are shear- 
periodic. Splitting off the background shear comes at the 
expense of an extra source term in the energy equation, 
which represents the work done by Reynolds stress due to 
the background shear (first term on the right-hand side of 
equation (Ig] ), see also [Stone &: Gardiner|[2010| ). 



3 NUMERICAL METHOD 

Since we are dealing with numerical convergence, we go into 
some detail explaining the numerical method used. We use 



operator splitting to solve the different parts of equations 
(13 1 - (16 1. Since we are interested in a balance between 



cooling and shock heating, it makes sense to use a Riemann 
solver to deal with the hydrodynamics. Throughout, we will 
use a uniform Cartesian grid covering the shearing sheet of 
size Lx X Ly, with grid spacings Aa; and Ay. 



3.1 Time step 

For stability, the time step is limited by the Courant- 
Friedrichs-Lewy (CFL) condition ( Courant et al.|1928 l: 



At — Co min 



Aa; 



Ay 



\u\ + Cs \V\ + Cs 



(17) 



where u = m/T., v = n/E, and the minimum is taken over 
all grid cells and Co is the Courant number. We have used 
Co = 0.4 throughout. Note that it is the perturbed velocity 
that enters the time step calculation. Because the large back- 
ground shear velocity has been split off, much larger time 
steps are possible ( Masset|2000a l. In low resolution studies, 
care must be taken that the time step does not become so 



large that two neighbouring rows are sheared apart ( Masset 
2000b I. This is never an issue for the resolutions adopted in 



this paper. 



3.2 Orbital advection 



Given a time step At, the most straightforward integration 



step involves just the first two terms of equations ( 13 1 - 
( 16|, which is linear advection at the background Keplerian 



shear ( |Masset|2000a|b| . In the following, we adopt the usual 
convention that X"j denotes quantity X evaluated at time 
index n and space indices i (in the x direction) and j (in 
the y direction). For each row at Xi, we therefore need to 
shift the solution by Sy = qi}xiAt. This can be achieved 
for quantity X by first shifting the solution by an integer 
number of grid cells A'^, where A'' is the nearest integer to 
5y/Ay. The rest of the shift Sy — NAy can be done to second 
order accuracy by standard methods (e.g. LeVeque]|2002 \ : 



-a{l-\a\){Gi 



(18) 



where a — 5y / Ay~N , a denotes a limited slope, and u = j — 
a/\a\ is the upwind direction. The second term on the right- 
hand side gives a first-order update, while the third term 
provides second-order corrections. Note that since N is the 
nearest integer to 5y/Ay, we always have that \a\ < 1/2, so 
that the above scheme is stable for any value of At. For the 
slope limiter, we use the same form as used in the Riemann 



solver (superbee, see Sect. 3.5.11 



Stone & Gardiner (20101 include the extra energy 



source term (first term on the right-hand side in equation 
(16l) in the orbital advection step. We found this can oc- 
casionally lead to negative pressures, and that stability is 
improved by instead integrating this source term during the 
X integration (see below). Furthermore, it proved numeri- 
cally advantageous to shift the pressure according to equa- 
tion ( |18| l rather than the total energy, again to avoid un- 
physical states when e is dominated by kinetic energy. 
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3.3 Cooling 

It is numerically convenient to split ofT the non-geometrical 



source terms and integrate them separately (Eulderink & 
Mellema 19951. In our case, this involves only the cooling 



source term, so in this step we solve 
de _ pfl 

dt'~{j-i)/3' 

or, since all other state variables (E, m and n) are constant 
during this step: 

dp pQ 



(19) 



(20) 



For the cases we are interested in. At <^ /3/Q, and we have 
found that a first order integration of equation ( |20[ | gave 
good enough results. 

3.4 Gravitational potential 

To calculate the gravitational potential $, we basically fol- 



low Gammie (2001 1. We first shift back the density in y to 
the time it was last periodic {t = tp). For this, we use the 
same algorithm as used for orbital advection. The resulting 
surface density is completely periodic in x and y, and the 
potential can then be found from 



$ = -2nGj2 ^ exp(ik • x - |k|5). 



(21) 



where Ek are the Fourier components of the shifted E and 
k = (fcx -|- qfl{t — tp)ky, kyY^ . This calculation can be done 
most effectively by using the Fast Fourier Transform (FFT). 
First, the Fourier components of the surface density are cal- 
culated using the FFT, the result of which is divided by 
|k|. An inverse FFT then gives us This potential is then 
shifted forward in y to the present time, after which the 
forces are obtained by taking finite differences. Even when 
using 5 = 0, the use of finite differences naturally intro- 
duces a smoothing length comparable to the grid scale for 
the gravitational force. 



3.5 Integrating in the x direction 

Taking only terms involving derivatives with respect to t and 



X from equations (13 1 - (16 1, together with the appropriate 



geometrical source terms ( [Eulderink fc Mellema 19951, we 
get: 



9E ^ dm 
dt dx 



dm 
~dt 



+ 



d_ 
dx 



0, 



m 



+ p] ^ 2Q.n - E 



dx ' 



dt 



d fmn\ , „, ^ 



(22) 



(23) 



(24) 



or the Reynolds stress due to the background shear (see 



also equation (12l). Note also that this would not hold if 



we had integrated the energy source term associated with 
the background shear during the orbital advection step as 



Stone & Gardiner (20101 



3.5.1 Roe solver 

Ignoring the source terms for the moment (we come back to 
these in section 3.5.2 below), the left-hand sides of the above 
equations have exactly the structure of ordinary Cartesian 
hydrodynamics, even though we are working with fluctua- 
tions on top of the background shear. Therefore, we can use 
standard techniques to deal with these. As mentioned above, 
since we are interested in shock heating, it makes sense to 
use a Riemann solver. This is advantageous not only when 
dealing with shocks, but for waves in general, because of the 
use of a characteristic decomposition of fluctuating quanti- 
ties. 

The above equations (still ignoring source terms) can 
be written concisely as dW/dt + dF/dx — 0, where W — 
(E, m, n, e)"^ and F = (m, m^/E -I- p, mn/T,, (e + p)m/Ti)'^ . 
Recalling that u = m/E, the Jacobian matrix A — dF/dW 
has eigenvalues u + c^, u — Cg and u, and eigenvectors 

\T 



ei = 

e2 = 

63 = 

e4 = 



(1, U + Cs,V,h + CsU) 



{l,u 
(0,0, 1,?; 



Cs,v,h- 

T 



CsU) 



u 
~2 



(l,u,« 

where w = n/E and 

7 



2 2 
U LU 

h= \ h - 

2 2 7 



1 E 



(26) 
(27) 
(28) 

(29) 



(30) 



A vector A — (Ae, Am, A„, Ae) can be projected 
onto the eigenvectors of A using the coefficients 



ai 



0.2 



04 



1 - 



1 



2cf 
Ae- 



+ 



7 - 



1 



2cf 

Ae - UA; 

A„ - vAf, 
^ \{h 



+ 



2 

dA„ 
2 



A„-f 



+ 



- wA„ 



2c, 



A„+ 



Am — uAp 



7 



A, 



Ae -I- uAm + vArA 



(31) 



(32) 
(33) 
(34) 



A linear hyperbolic system of the form dW/dt + 
AdW /dx — can be solved by projecting the state differ- 
ence between neighbouring grid cells onto the eigenvectors 
of A, or, in other words, decomposing the state difference 
into characteristic waves: 



OkSk, 



(35) 



de d 
dt dx 



, „ mn d<^ 

(e-fp)— ) = qVL— -m- 



e; -e -dx- (25) 

Note that all source terms in the energy equation are due 
to the kinetic part of the total energy: the thermal energy 
does not change directly due to work done by self-gravity 



where we have omitted the subscript j for clarity. The state 
at the cell interface is found by subtracting all right-going 
waves from W^: 



-1/2 



Wj - ^ max(sign(Afe), 0)afeefc, 



(36) 
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where Afc is the kih eigenvalue of A. The interface flux 
F,_i/2 = AW,^i/2: 



sign(Afc) + 1 



akXkSk- 



(37) 



Alternatively, the interface state can be found by taking into 
account all left-going waves from Wi_i: 

Wi_i/2 = W,_i - ^ min(sign(Afc), 0)afeefc, (38) 

k 

from which we can again obtain an interface flux 

-akXkSk- (39) 



•P w sign(Afc) - 1 

* «-i/2 - * i-1 - Y 



Combining both expressions for Fi_i/2, we get 
F,_i/2 = ^ ^Fi_i + Fi - ^ |Afc|afcefe j , 



A flrst-order state update is then given by 

At „ X 

Aa; 



-1/2 



-1/2) 



(40) 



(41) 



Equations ( |22[ |-(|25[l are in fact non-linear. However, 
Roe ( 1981) introduced a suitable linearisation, in which the 
matrix A is replaced by a matrix that is averaged over two 
neighbouring grid cells. By demanding that the resulting 
solution should be exact if the two cells were connected by 
a single shock wave, |Roe| ( |1981[ ) found that the eigenvec- 
tors (26l-(29l and projection coefficients (31|-(34| should 



be evaluated at the so-called Roe-averaged state 



-(- \/Ei_i 
and an average sound speed 



h = 



= (7- 1) ( ^ - - 



(42) 
(43) 
(44) 

(45) 



We can then use the flux function (401 to update the state 



using ( 41 1 



One can obtain a method that is second-order accurate 
in both space and time by using modified projection coeffi- 
cients puldcrink fc Mellema|1995[ [LeVcquc, 2002j ) 

At, 



fflfc 1 



Ax 



(46) 



However, near discontinuities any method that is more than 
first-order accurate will introduce unphysical oscillations 
( Godunov|1954 1. It is therefore necessary to introduce a flux 
limiter (h: 



a.k = ak 



irk)^JM 



(47) 



where rk ~ ak/aku is a parameter comparing Ofe to the value 
of flfe in the upwind direction. In smooth flow, we expect 
rfc ~ 1 and we obtain a second-order method if 0(1) = 1. 
The functional form of S is chosen so that no oscillations 



are introduced by the term proportional to 0, which can 
be made mathematically precise using the concept of To- 
tal Variation (TV) (e.g. LeVcquo 2002). Any chosen limiter 
function must by TV Diminishing (TVD) in order for it not 
to introduce oscillations near shocks. A fairly general class 
of TVD limiters can we written as 



— max(0, min(l, s9),mm{s, 0)), 



(48) 



which is TVD for 1 ^ s ^ 2. The case s = 1 corresponds to 
the minmod limiter, which is the most diffusive choice, while 
s = 2 corresponds to the superbee limiter, which in general 
gives the sharpest shocks. Because of the use of a flux lim- 
iter, no explicit artificial viscosity is needed to stabilise the 
numerical scheme. This however comes at the price of giving 
up second-order accuracy in non-smooth regions of the flow. 
In all simulations presented, we have used the superbee flux 
limiter. 



3.5.2 Source terms 

The system dW /dt -f dF/dx = S can be solved by com- 
bining the Riemann solver solution to dW /dt -\- dF/dx = 
with either a solution to dW /dt = S or dF/dx = S. The 
latter choice, called stationary extrapolation (Eulderink & 



Mellema 19951, is the preferred choice when a balance be- 



tween source terms and flux gradients is expected to arise 
in the solution. This is the case for example in global simu- 
lations of Keplerian discs ( Paardekooper &: Mellema 2006| , 
where in the radial direction there exists a balance between 
gravity, the centrifugal force and pressure. In the present 
case, we are dealing directly with perturbations. It is then 
more advantageous to use dW /dt — S. 

In particular, it is possible to deal with epicyclic oscil- 
lations in a way as to conserve the energy associated with 
epicyclic motion to round-off error dGardiner fc Stone|2005| 



Gressel fc Ziegler|2007| [Stone fc Gardiner|2010| ). Recall that 
we are solving 



dm 
dt 
dn 



{q-2)nm, 



(49) 
(50) 



dt 

where Sx ~ —"Sd^/dx is the source term due to self-gravity, 
together with d'E/dt — dp/dt — 0. The energy associated 
with epicyclic motions can be conserved by integrating these 



equations using a Crank-Nicholson scheme (Stone & Gar- 
diner|2010| |: 



2Atnn -\- AtSx, 
{q - 2) AtQm, 



(51) 
(52) 



where x — (x" -\- x"')/2 is a time average, and [x] is short 



for X 



n + l 



x". Solving these two equations gives: 



2Atn (2n" - AtQ. (2 -q)m" + S^/^) 



-2Atn (2 - q) 



2 -\- Af^Q.-^ (2 - g) 

m" + AtSx/2 + Atfln" 



(53) 
(54) 



2 + Af^Q? (2 - g) ' 

Note that Sx = ^(E), which is therefore constant during 
the source term integration. 

As an alternative, Sx could be left out in this step, 



and integrated using stationary extrapolation ( Eulderink & 



6 S.-J. Paardekooper 



Mellema 19951. This could be advantageous when objects 
form in which pressure is balanced by self-gravity. This bal- 
ance would then be recognised by the Roe solver, resulting 
in no evolution away from this steady state. 



3.5.3 A full time step 

A complete integration step for the x direction then consists 
of the following steps: 

• Prediction: evolve the momenta under influence of the 



source terms using equations ( 49 1 and ( 50 1 for half a time 
step. A simple Euler integration is sufficient in this step. 

• Riemann solver: use the output of the previous step to 
calculate a state update dWhydro for a full time step. 

• Source integration: use equations ( |53[ ) and ( |54[ ) calcu- 
late the source update to the state dWsourco, keeping E and 
p constant. For the state variables on the right-hand side we 
use W" -I- dWhydro/2. Note that we need & for both E" 
and E"+^, so that we need to recalculate the gravitational 
potential after the previous step. 

• Update the state according to W"+^ = W" -|- 
dW 

source ro . 



3.6 Integrating in the y direction 

For the y-direction, we proceed in a similar way. The gov- 
erning equations in terms of m and n read: 



9E 9n _Q 
dt dy 



dm d /'mn\ 
-Ot+0-y{^)='' 



+ p] =-E 



dn d f n? 



de 9 / , > n 

9t + a^((^ + ^)E 



3.6.1 Riemann solver 



9$ 

dy ' 



'■IT- 
dy 



(55) 



(56) 



(57) 



(58) 



The left-hand sides of above equations have the same struc- 
ture as their x-equivalents, which leads to similar expressions 
for the eigenvectors and projection coefficients: 



ei = {l,u,v + Cs,h + Cbv) 
e2 — {l,u,v ~ Cs,h — Csv) 
eg = {0,l,0,uf, 

2 2 

e4 = {l,u,v, Y + y) > 



(59) 
(60) 
(61) 

(62) 



and 



ai = 



7^1 
2c2 



a2 = 



7-1 
2c| 

A, - 



as = A„ 

7 - 

a4 — — 





+ y 


uAm 


- vAr 




+ y 


uAm 


- vAr 


uAp 






2 

U — 



A„+ 



An — vAp 



A„+ 



An — vAp 



h — u — V ^ Ap — Ae + uAm -f uA„] 



(63) 



(64) 
(65) 
(66) 



Again, using Roe-averaging and flux-limiting, these formulae 
can be used to construct a second-order correct update for 
the state vector. 



3.6.2 Source term integration 

There are no epicyclic oscillations when considering the y- 
direction only, and we can therefore adopt a simple integra- 
tion scheme 



n + l 



= n" -f AtSy, 



(67) 



where 5*1, = — E9"l?/9y, together with dE/dt = dm/dt = 
dp/dt — 0. Alternatively, stationary extrapolation could be 
used for this step. The complete integration step is equiva- 



lent to that in the a;-direction (see section 3.5.3 1 



3.7 Backup fluxes 

We expect strong density and pressure contrasts to arise in 
the simulations, especially in cases that fragmentation will 
happen. These conditions provide a challenge for numerical 
methods, and it is important to assess where the method 
could fail. Since Riemann solvers require the use of the to- 
tal energy, they are more likely to fail (i.e. predict negative 
pressures) when |uj ^ Cs ( Einfeldt et al.|1991 l. 

The maximum possible time step restricted by the CFL 
condition is QAt ~ QAx/cs — Ax/H. Provided that H is 
properly resolved (as it should be), we have that Q,t <^ 1. 
For P ^ 1, we therefore expect the simple Euler integration 
of the cooling not to give unphysical states. 

The source term integrations during the x and y inte- 
gration steps all involve dT,/dt = dp/dt = 0. This means 
that the states fed into the Riemann solver by the predic- 
tion step (see section 3.5.3 1 are all physical (i.e. E > and 



p > Q). The only place where the scheme can break down is 
therefore the Riemann solver itself giving negative pressures 
or surface densities. This can for example happen when the 
Riemann problem between two cells is not linearisable ( |Ein-| 
feidTet al. 1991]. 



A possible way of dealing with such failures is to add 
more diffusion by returning to the first-order fluxes locally. 
Although the first order Roe fluxes are readily available, 
they can still produce unphysical states if the Riemann prob- 
lem can not be linearised. As an alternative, we use the 



et al. 



1983 



feldt et al. 



more diffusive, Harten-Lax-van Leer (HLL) solver (Harten 
. This scheme is in fact positive definite (Ein- 
199T| ), which means that, given physical input 



states, it will never lead to negative pressures or densities. 
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4.1 Epicyclic motion 

In the absence of self-gravity and any gradients in x or y, 
there exist oscillating solutions to the governing equations: 



, , 2r2no . 
m(t) — mo cos Kt H sm Kt, 

K 

n(t) — nn cos Kt — sm Kt, 



(72) 
(73) 



where k^ = 2(2 — q)i^^ is the square of the epicyclic fre- 
quency. In a Keplerian disc, with q — 3/2, we have that 
K = Q. 

We take a grid of 128^ -1/2 s: {x,y) ^ 1/2, take 
Cs — 0.01 and give the whole grid a velocity perturbation 
of 11 = O.lCs. The resulting evolution of u is shown in Fig. 
[l] Because of the Crank-Nicholson time integration, the am- 
plitude remains constant up to round-off error. The phase 



accuracy is determined by the time step ( Stone & Gardiner 
2010[ ). For the case of Fig.[l] flAt = 0.13, which leads to a 



phase error of less than 1% over the course of the simulation. 



Figure 1. Epicyclic motion in tlie absence of any gradients on a 
grid of 128^, -1/2 ^ {x,y) (i 1/2, and Cs = 0.01. 



The corresponding flux function is, in the x-direction, given 

by 

F,-i/2 = b+-b- ' ^ > 



where 



b = min(0, Afc, Afe,j_i), 

k 

fo"*" = max(0, Afe, Afc,i) 

k 



(69) 
(70) 



are measures of the minimum and maximum possible wave 
speeds encountered in the Riemann problem. Here, Xk,i are 



the eigenvalues of the Jacobian matrix (see Sect. 3.5.1 1 based 
on the state and flux of cell i, and A^ are the eigenvalues 
based on the Roe-averaged state between cells i and i — 1. If 
the Roe flux is found to lead to an unphysical state, the flux 
is replaced by the HLL flux. Usually, this is only necessary 
once in every 10^ updates. 



3.8 Boundary conditions 

Periodic boundary conditions are used in the y direction. In 



4.2 Linear shearing waves 

A more challenging problem, that also tests the self-gravity 



solver, consists of evolving a linear shearing wave ( Gammie 



20011. Below, we first derive the governing equation, basi- 



cally following Gammie ( 1996 1, but adapted to our notation 



and neglecting viscosity and magnetic fields. 



4-2.1 Governing equation 

Consider linear, adiabatic perturbations of the governing 
equations, so that the pressure perturbation pi — CgEi, 
where Cg is the adiabatic sound speed. The continuity equa- 
tion and the two momentum equations then read 







1 DEi 
Eo Dt 


dui 


^ dvi 






dx 


dy 




Dui 


- 2Q.VI + -|- 
Eo 


dpi 






iDt 


dx 


dx 




+ (2 


- g)r2iti -1- ^ 


dpi 


dy 


Ijt 




dy 



0, 



= 0, 



= 0, 



(74) 
(75) 
(76) 



where D/ Dt = d/dt — q^lxd/dy is the convective derivative 
with respect to the unperturbed flow. Note that D/Dt and 
d/dy commute, but that 



the 


X direction, the sheet is shear-periodic; for example at 






.A 


fDf\ 


the 


inner boundary we have that 




idi) 


dx 


\ Dt) 



+ qVl 



dl 
dy' 



(77) 



E(a;,y) = E(a; + L.j:,y - qQ,Lxt), 



(71) 



where Lx is the size of the sheet in the x direction. The 
required shift in y is again performed using the same method 



Take the convective derivative of equation (741 
1 D^Ei 



as used for orbital advection (see Sect. 3.2 1 



1 D^Ei 



+ 



Eo Df^ 
d ( Dui 



D_ 
Dt 



dui 
dx 



+ 



Eo Dt^ dx \ Dt 



dvi 
dy 

d_ fDvi\ 



D_ 
Dt 



(78) 



4 TEST PROBLEMS 

The Riemann solvers were tested using standard one- 
dimensional shock tubes and two-dimensional Riemann 
problems. Below, we describe two test problems that are 
speciflc to the shearing sheet. 



and insert equations ( 75 I and ( 76 1 
1 D^Ei , d 



Eo Df^ ~^ dx 

d 



2nvi 



dy 



2{q - l)nui 



1 dpi d^i 
Eo dx dx 
1 dpi d^i 
Eo dy dy 



+ 



= 0. 



(79) 
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The potential vorticity ^ is given by 
{2-q)n + dv/dx-du/dy 



(80) 



(81) 



which can be approximated by a background value 
(2 — g)Sl/Eo and a perturbation 

^ _ dvi/dx-dui/dy _ {2-q)n Si 
Eo So Eo 

For our test problem, we are interested in adiabatic pertur- 
bations, for which potential vorticity is conserved. Within 
the linear approximation, this means that ^ is a shearing 
wave with constant amplitude, or D^/Dt = 0. We must 
therefore have that 



dui dvi El 

-7r~ ^ ~^ 2 - q)n- Eo^i, 

ay ox Eo 

with ^1 a combination of shearing waves: 
Ci = ^ {x) exTp(iqQkytx + ikyy). 



(82) 
(83) 
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We can use equation ( 82 1 in equation ( 79 1 
d ( 

2qnvi 



1 D^Ei 



+ 



1 dpi 9$i 
Eo dx dx 

2{q ^ 1){2 ~ q)n^ 



Eg Dt'^ dx V 
1 d'^pi _ 

Eo dy^ dy^ Eo 

-2(5- l)nEoCi = 0. (84) 

Differentiating the continuity equation with respect to y, 
and using equation (821, we obtain an expression for vi: 

a^^ + ay^j"^^(2~g)^^^^ 

At this point, we decompose the solution into shearing 
waves: 



Xi = Xi{t) exp (i(k^{t)x + kyy)) , 



(86) 



with kx{t) = kxfi+q^kyt. For notational convenience, we set 
Xi = Xi, the exponential factor being taken as read. Note 



that now ^i is a constant. From equation (85 1 we obtain 

(87) 



2 iky dEi . ^o^l -7 ^ f 

k VI = — — ikx[2 - q)il— ikxT.0^1, 

Zjo (it Zjq 



with k — kx{t) + ky, while equation (84 1 reads 

^ + 2qnik^vi + k% + e^i- 

Eo dt'^ Eo 

2{q - 1)(2 - g)n^^ - 2{q - l)fiEo6 = 0. 
Eo 

Using the expression for vi above, pi = CgEi, and <l?i 
— 27rGEi/fc, we finally obtain 

1 d^El ^ ^k,,ky 1 dEi , 

— 2qQ, — — + 



(88) 



Eo df^ 



fe2 Eo dt 



U2\ y 

2(2 - q)n^ + k^cl - 2TTGEok - 2q{2 - qp^ -i + 

k-' J Eo 



2n{l-q-^] EoCi = 0. (89) 



This equation is equivalent to equation (15a) of |Gammie| 
(19961, but without viscosity and magnetic fields. 



Figure 2. Evolution of a non-selfgravitating linear sliearing wave. 
The tliick solid curve indicates the solution to equation (jsojl , while 
the other curves indicate results from hydrodynamic simulations 
at different resolutions. 



4-2.2 Numerical results 
r a shearing \ 

with Eq = 1/40 and Q = 1 of initial amplitude Ei/Eo — 
0.0005 with k^ = -Atv and ky = 2n. |Gammie| (^2001 ) con- 
sidered a similar problem, but with pi = initially, so that 
equation ( |89[ l does not apply since the initial perturbation 
is not adiabatic. 

In Fig. [2] we consider the case without self-gravity for 
different numerical resolutions. A fundamental scale to re- 
solve is the scale height H = Cs/Jl ~ 0.07. At the lowest 
resolution, H is resolved by two grid cells only, which leads 
to strong diffusion of the wave. Since Riemann solvers ac- 
tively use the sound speed, or, more general, the physical 
scales in the problem, to compute fluxes, not resolving the 
physical scale of the problem can lead to more excessive 
diffusion than for other numerical methods at comparable 
resolution. For eight grid cells per scale height (A'^ = 128), 
good agreement with linear theory is obtained up to Qt = 8. 

We now turn to the case with self-gravity in Fig. [3] us- 
ing 5 = 0. For N = 128, the numerical result follows linear 
theory in a similar way as for the case without self-gravity, 
which also compares well with (Gam mie||2001| his Fig. 1). 
Again, not resolving the physical scales of the problem leads 
to excessive diffusion for A'^ — 32, apparently more so than 
m |Gammie| ( [200T| ). This can be attributed to the fact that a 
Riemann solver actively uses the physical scale when com- 
puting the fluxes, so a penalty is paid if the necessary scales 
are under-resolved. 



5 INITIAL CONDITIONS 

Initial conditions require special attention. In order for a 
steady gravito-turbulent state to be set up, we have to avoid 
initial transients ( Paardekooper et aLpOll l. Even though in 
local simulations no edges can form due to a radial depen- 
dence of the cooling time scale, it still takes a finite time 
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Figure 3. Evolution of a self-gravitating linear shearing wave. 
The thick solid curve indicates the solution to equation \S9\ , while 
the other curves indicate results from hydrodynamic simulations 
at different resolutions. 



for the turbulence to develop (fifdcvci ~ 10). For cooling 
times smaller than or comparable to tdcvoi, the disc will be 
cooled down to Q < 1 before the disc can start to balance 
the cooling, which can trigger fragmentation even in cases 
where a steady gravito-turbulent state may exist. We there- 
fore choose to hold Q > 1 until fit = 50. This is long enough 
for turbulence to develop and keep the disc from fragment- 
ing artificially. After Qt = 50, the disc is free to cool down 
and fragment. 



We take Sp = 1/320 with 



similar to Gammie (20011. This makes H 



The highest resolution considered by Gammie (20011 was 



1 and (5 = 1, 

= cs/n » 0.01. 



Ny 



1024, resolving H by approximately 10 grid 



cells. This will be our standard resolution, and we go up by 



factors of 2 from there. To compare with Gammie (20011, 
we use 5 = unless otherwise specified. The initial velocity 
field is seeded with white subsonic noise to let the turbu- 



lence develop. Following Gammie (2001 1, we take 7 = 2. We 



consider 16 different cooling times, varying (3 between 1 and 
50. Each simulation is run until Qt = 1000. Since we will 
find that fragmentation can be a stochastic process, we run 
four versions of a simulation, where we vary the phase of the 
initial seed noise, keeping the amplitude constant. 



6 RESULTS 

The aim is to test for numerical convergence of the deter- 
mination of the critical cooling time scale /3c. We do this by 



running simulations similar to those in Gammie (2001 1, but 



at higher resolution and for longer time spans. 



6.1 Reproducing previous results 

First of all, we try to reproduce previous results at the stan- 
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Figure 4. Measured a parameter as a function of imposed cooling 
(open circles) together with the prediction of equation ([sjl (solid 
line), for = Ny = 1024 and Htmax = 100. The vertical dotted 
line shows the fragmentation boundary. 



usually done, we define the disc to have fragmented when an 
overdensity of IOOSq survives for several cooling time scales 
(e.g. |Meru fc Batel[20TT] |Rice et al.|2"oTT| ) . 

We calculate the total stress in the sheet in the usual 



way (e.g. Gammie 2001 Rice et al. 20111. The average 



Reynolds stress is given by 



(90) 



where () denotes an average over the whole computational 
domain. The gravitational stress is most easily determined 
in the Fourier domain: 



{Gxy} — ^ 



TvGk^ky |Eki 



(91) 



The total stress can be parametrised using the a- 
prescription: 



(Hxy) + (Gxy) 



(92) 



which can then be compared to equation ([3|. We average the 
measured values of a over QAt = 20 to get a single value 
for a given simulation. 

In Fig. |4] we show the measured values of a together 
with the prediction of equation (|3|. Over the range of /3 we 
consider, we find agreement to within 5%. Simulations with 
(3 < 4 were found to fragment, which is in good agreement 



with Gammie ( 2001 1, who found (3^ = 3. Moreover, the max- 
imum value of a the disc can sustain is Omax = 0.054 (for 



(3 = 4) , in good agreement with Rice et al. ( 2005 1 , who found 
Qmax ~ 0.06. We also note that the measured rms density 



fiuctuations agree with the results of Cossins et al. (20091 



dard resolution (Nx 



Ny 



1024) and for Qt < 100. As is 



Since the physical scale of the instability, the most un- 
stable wavelength At ~ is resolved, one might argue that 
these results should be converged with respect to numerical 
resolution. Moreover, Gammie (2001) showed that the mea- 
sured value of a is independent of resolution if A*' ^ 512. 
We have confirmed that the results depicted in Fig. |4] do 
not change when decreasing the resolution by a factor of 2. 
However, the value of a is not necessarily a good indicator 
of numerical convergence. Given the prescribed amount of 
cooling, the disc will try to generate enough heating to make 
up for the energy that is removed. In the present set-up, it 
can only do that by generating the necessary stresses. Un- 
less the simulation is dominated by numerical viscosity, the 
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Figure 5. Evolution of ttie maximum surface density (top panel) 
and the measured value of a (bottom panel) for a run with A^^ = 
Ny = 1024 and 13 = 5. The dotted li 
indicates the prediction of equation 



measured value of a will always be very close to the predic- 
tion of equation ([3|; otherwise, the disc can not maintain a 
steady state. This, however, does not necessarily mean that 
the result makes sense, physically. In particular, convergence 
with respect to a does not imply convergence for the value 
of /3c. 

In a similar way, resolving the physical scale of gravita- 
tional instabilities is only a necessary condition for numerical 
convergence. It would probably be sufficient if no other pro- 
cesses were going on, and if the evolution is predominantly 
on dynamical time scales. If very slow time scales on small 
scales are involved, it is likely that higher resolutions are 
required to capture the numerical evolution correctly. We 
will see below that processes happening on the cooling time 
scale are critical in determining whether the disc will frag- 
ment or not. It is therefore expected that for increasing /?, 
higher resolutions are required. 



6.2 Longer time spans 

We now keep the resolution fixed at A^^, — Ny — 1024, but 
integrate all simulations until Jltmax ~ 1000 (or until frag- 
mentation occurs). In Figs. [s] and [o] we focus on the case of 
/3 = 5, which was found not to fragment for Qt < 100. 

In the top panel of Fig. [5j the evolution of the maxi- 
mum surface density in the sheet is shown. Before fit — 250, 
Smax/So stays below 50, approximately, and the measured 
value of a agrees with equation (|3| (bottom panel of Fig. 
[5|. The top panel of Fig. |6] shows a snapshot of the surface 
density at Qt — 100. This looks like a very good example 
of steady gravito-turbulence, with density fluctuations that 
are consistent with those found by Cossins et al. (20091. 

However, after fit = 250, something interesting hap- 
pens. Suddenly, the maximum surface density shoots up 
to values above 100, indicating fragmentation. The bottom 
panel of Fig. [6] shows a snapshot of the surface density at 
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Figure 6. Surface density, in terms of the initial surface density, 
on a logarithmic scale, for a simulation with Nx = Ny = 1024 
and /3 = 5. Top panel: fit = 100, bottom panel: fit = 400. 



Qt — 400, after the disc has fragmented. Only a single frag- 
ment was formed around Qt = 250, in contrast to simula- 
tions with P < Pc which usually show ~ 5 — 10 fragments, 
initially, which can subsequently merge. 

The reason for this fragmentation at high values of /3 
lies in the nature of the gravito-turbulent state. Even before 
true fragmentation occurs, clumps are formed and destroyed 
on a continuous basis. This can be appreciated from the top 
panel of Fig. [5j where the peaks in Emax indicate a clump 
being destroyed. The root-mean-square density fluctuation 
is of order unity, while the maximum surface density reaches 
values of Emax/So — 50 several times. One clump that does 
not make it to collapse can be spotted near x = 0.05 and 
y = —0.45 in the top panel of Fig. [6] 

Clumps of size ~ H can survive the tidal shear if their 
size is less than the size of their Hill sphere. If we take the 
surface density within the clump to be constant for simplic- 
ity, we must have that 



H <Ro 



3A/. 



1/3 



(93) 



where Rq is the radial distance to the central star and A/, is 
its mass. This condition can be recast in terms of the local 



Stochastic disc fragmentation 11 




400 600 
£2t 



Figure 7. Evolution of tiie maximum surface density (top panel) 
and the measured value of a (bottom panel) for four realisations 
with Nx = Ny = 2048 and /3 = 9. The dotted line in the bottom 
panel indicates the prediction of equation (jsjl. 



value of Q: 
1 

3' 



Q < 



(94) 



In other words, keeping the temperature fixed, we only need 
an increase in surface density of a factor of 3 over the back- 
ground Qo ~ 1. state to form a clump that can resist the 
shear. Once formed, these clumps will in general contract on 



a cooling time scale (Kratter & Murray-Clay 20111. Their 



survival depends mainly on if they can resist the weak shocks 
that sweep around in gravito-turbulence. Since shock heat- 
ing is very localised, this makes fragmentation a stochastic 
process: there will be a large spread in clump survival times, 
until the first lucky clump survives long enough for collapse 
to proceed. It should be noted that the condition given by 



equation ( 93 1 is not necessary if the cooling time scale is 
comparable to the dynamical time scale. If cooling acts on 
a dynamical time scale, there is no time for the clump to 
shear apart before it collapses. 

We have observed fragmentation up to /3 = 7, more 
than twice the critical cooling time scale found by |Gammie| 
1 2001 1. The corresponding maximum value of the stress is 
Omax ~ 0.03. For larger values of /3, the disc remained in a 
steady, gravito-turbulent state for Qt < 1000, with values of 
a that agree well with equation ([3|. 



6.3 Higher resolution 

We find that increasing the resolution by a factor of 2 {Nx ~ 
Ny — 2048) leads to easier fragmentation at higher values of 
13. As an example, we show in Fig.[7]four simulations at 13 — 
9, differing only in the phase (not magnitude) of the initial 
noise. Two of the discs fragment, one at Qt « 500 and one 
at Qt ~ 750. The other two discs maintain a steady gravito- 
turbulent state for the full length of the simulation. This 
nicely illustrates the stochastic nature of disc fragmentation 
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Figure 8. Surface density, in terms of the initial surface density, 
on a logarithmic scale, for simulations with Nx = Ny = 2048. 
Top panel: /3 = 20, fit = 160, bottom panel: /3 = 40, Ut = 370. 



at high values of /3: only in two out of four simulations does 
a clump survive for long enough for collapse to proceed. It is 
expected that if the simulations would be continued, in the 
end all of the four realisations should show fragmentation. 
Note that fragmentation now occurs for three times higher 
values of [3 than in Sect. |6.l| 

We have found clumps that can survive shear to form 
for all values of /3 we have considered (/3 5C 50). In Fig.js] two 
examples are shown of clumps with E/Eo ~ 10 for (3 = 20 
(top panel) and j3 — 40 (bottom panel). However, while 
clumps form readily in all simulations, the vast majority do 
not survive. In the top panel of Fig. [7| at least 10 clumps 
were formed and destroyed before the first disc fragments. 
In discs that do not fragment, more than 20 clumps form 
during the time span of the simulation, but none of them 
collapse into bound fragments. Given this low success rate, 
simulations should span at least 100 cooling time scales to 
capture these events. This becomes impractical for very high 
values of (3. The highest value of (3 for which we have found 
fragmentation is /3 = 20 (with Nx ~ Ny — 4096), almost 7 
times higher than what was found in Sect. |6.1[ 

The reason why this behaviour can not be captured at 
low resolution, lies in the fact that we need the clumps to 
survive for QAt ~ /3. It is not enough just to resolve the 
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Figure 9. Evolution of the maximum surface density for simula- 
tions with /3 = 8 in a domain that is reduced in size by a factor 
of 8, for two different resolutions. 



length scale H, the numerical scheme needs to be able to 
maintain a coherent clump of size H over many dynamical 
time scales. The resolution required will depend on the de- 
tails of the numerical implementation, and will increase for 
larger values of /3. In addition, long time integrations are 
needed to weed through all the clumps that fail to collapse. 

The need for high resolution can be appreciated fur- 
ther by looking at smaller domains, so that only a single 
transient clump is present at any time. The evolution of the 
maximum surface density is then linked directly to the evo- 
lution of this clump. In Fig. [9] we show the evolution of the 
maximum surface density for /3 = 8 in a domain that is re- 
duced by a factor of 8 compared to the standard run (i.e. 
L^^ Ly ^ 1/8). Note that still H L^, and that iV = 128 
now corresponds to our standard resolution. Neither of the 
two simulations shows fragmentation, but transient clumps 
can be clearly identified. In the high-resolution run, it takes 
approximately ilt — 40 for a clump to reach its maximum 
density, which corresponds to 5 cooling time scales. In the 
low resolution run, clumps do not survive this long, leading 
to less-pronounced surface density peaks. Therefore, clump 
survival appears to be linked to numerical resolution. 

The frequency at which clumps appear decreases for 
higher values of /?. While for P — 9, approximately 20 failed 
clumps can be identified over a time span of fit — 1000, 
for /3 = 50 only ~ 4 failed clumps appear. The frequency 
appears to be roughly proportional to 1//3. Together with 
the fact that the clumps contract on a cooling time scale, 
this makes fragmentation very rare, but not impossible, for 
higher values of /3. 



6.4 Additional numerical effects 

In this section, we discuss the infiuence of two more numer- 
ical parameters: the smoothing length 5 and the flux limiter 
(j). Previous studies | |Gammie]|2001| [Rice et al.||201l| used 
S — 0, which, at high resolution, allows for variations in the 
potential on scales much smaller than a scale height H. If 
we interpret the two-dimensional approximation as resulting 
from vertically averaging the three-dimensional equations, 
we do not expect to find such small scales. If these scales 
are important, this calls for fully three-dimensional simula- 
tions. We have performed additional runs using 5 = Cs/f2, 
where the average sound speed was calculated every time 
step. Note that this is a fairly large value for 5; in two- 
dimensional disc-planet interaction studies, the gravitational 




Figure 10. Evolution of the maximum surface density (top 
panel) and the measured value of a (bottom panel) for four reali- 
sations with Nx = Ny = 2048 and /3 = 7 with 5 = H. The dotted 
line in the bottom panel indicates the prediction of equation ([3| . 



potential of the planet is usually smoothed over a distance 
H/2 (e.g. |Baruteau fc Masset|2008l ). Moreover, the true disc 
thickness will be smaller than Cs/fi because of self-gravity. 
Using S = Cs/O may therefore introduce more smoothing 
than necessary. 

The results for /3 = 7 are shown in Fig. |10| As in the 
case with 5 = 0, fragmentation appears to be stochastic, 
with one out of four realisations forming a bound fragment 
within fit = 1000. This effect therefore does not rely heav- 
ily on small scales in the gravitational potential. However, 
no fragmentation was found for /3 > 10, in contrast to the 
simulations with 5 = 0. Thus, while the overall picture is 
very similar, with transient clumps clearly visible in Fig. 
|10| the efficiency of stochastic fragmentation is reduced for 
S = Cs/fl. Three-dimensional simulations are necessary to 
see which case is appropriate. 

A similar story holds for changing the flux limiter. Sim- 
ulations using the diffusive minmod flux limiter, obtained by 
setting s = 1 in equation (481, still show transient clumps 



and stochastic fragmentation, but only up to /? = 5 for the 
standard resolution. This is expected: since the minmod lim- 
iter introduces more numerical diffusion, it is more difficult 
for transient clumps to survive long enough to collapse. It 
is worth pointing out that in some multi-dimensional prob- 
lems, the minmod limiter appears to give more diffusion than 
finite difference codes ( [Paardekooper et al. 12008] ). 



7 DISCUSSION 

We have shown in the previous section that disc fragmen- 
tation is a stochastic process in the simplified system un- 
der consideration. It is likely that some of the simplifica- 
tions made will affect the efficiency of disc fragmentation 
for longer cooling times. We discuss the three most impor- 
tant ones below, before we focus on possible implications. 
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7.1 Limitations 

First of all, we have worked in two dimensions only. This 
means that all quantities should be thought of as being in- 
tegrated over the disc scale height H. While this is straight- 
forward for the gas density, velocities and pressure, the po- 
tential due to self-gravity poses a problem. The force due to 
self-gravity should be smoothed over a length comparable 
to the scale height. While introducing a smoothing length 
did not change the results in a qualitative way (see Sect. 
|6.4[ ), three-dimensional simulations are needed to determine 
whether the smoothed or the unsmoothed results are more 
appropriate for three-dimensional discs. 

Second, the use of the simple cooling law of equation 
([2| is questionable as soon as fragments form. When the 
density goes up by orders of magnitude, cooling should slow 
down dramatically. This can be a very important effect in 
the present situation, where clumps are formed on a continu- 
ous basis, trying to cool down before they are destroyed. The 
survival of clumps strongly depends on their ability to cool, 
and it is likely that when cooling slows down at increasing 
surface density, the likelihood of survival goes down. 

Finally, we have considered local models only. While 
this is advantageous for studying a steady gravito-turbulent 
state, it is impossible to say how the disc came into this 



state (e.g. Kratter & Murray-Clay 20111 or what happens 



to any fragments that have formed. For example, it is likely 
that these newly-formed objects are subject to rapid radial 
migration ( Baruteau et al.|2011 1. Also, any interaction with 
global modes ( [Lodato &: Rice||2005[ ) is necessarily excluded 
from the current study. 



7.2 Implications 

The most important conclusion from the results presented 
in Sect. |6] is that there is no such thing as a steady gravito- 
turbulent state in the simple case of /3-cooling, at least not 
for P < 20. Even though the disc can balance the imposed 
cooling for many dynamical time scales (see Fig.[7|, there is 
always a chance the disc will fragment at some point. This 
means there is no rock-solid criterion for disc fragmentation 
based on the cooling time scale. Fragmentation just gets less 
likely for higher values of p. The classical cooling criterion, 
which states that cooling should occur on a dynamical time 
scale, is in a way a condition for fragment survival as well: 
for /3 ~ 1, clumps can not be sheared apart by tidal forces 
(see Sect. 



2 I, ensuring essentially a 100 % survival rate. 



The maximum time span we have considered is fitmax ~ 
1000. Since the self-gravitating phase of protopl anetary discs 
is thought to last only for ~ 10^ years (e.g. Laughlin & 



Bodenheimer 19941, if we take our sheet to be located at 

max available before the disc is 



100 AU, there is less than Qtn 
no longer self-gravitating. It has to be kept in mind that even 
though steady gravito-turbulence may not exist formally in 
the simple case of /3-cooling, it is only necessary to be steady 
for a finite time span. There is therefore a value of /3 for 
which fragmentation becomes impractical, because it just 
takes too long for a rare event (the survival of a clump) to 
happen. It is difficult to say, for a given value of j3, exactly 
when the disc will fragment. This will for example depend on 
the total area of the sheet under consideration; larger sheets 
have a better chance to fragment. Since the vast majority of 



clumps will not survive, and since they contract on a cooling 
time scale, it is necessary to integrate for several lOO's of 
cooling time scales to see fragmentation. 

Based on these results for discs with simplified thermo- 
dynamics, it is dangerous to model the self-gravitating phase 
of disc evolution using an a-model. As soon as fragmentation 
sets in, stresses become dominated by the fragment (see for 
example the bottom panel of Fig. [5|. The extent to which 
the fragment can come to dominate will depend on its fi- 
nal mass, which is likely not to be captured very well in 
the current simulations because of the simple cooling law. If 
cooling becomes less efficient at higher density, the growth 
time of the fragments will likely go up, and the impact of 
the fragment will not be as dramatic as depicted in Fig. 
[5] Moreover, it is likely that the fragment will be subject 
to rapid orbital migration (Barute au et al.|2011 1, leaving its 
place of birth and perhaps allowing the disc to resettle into a 
gravito-turbulent state, albeit at lower mass. However, this 
all depends very sensitively on the mass evolution of the 
fragments, which is poorly constrained at the moment (see 
e.g. [Kratter et al.||2010| |Boley et al.|2010 |. 

Fragmentation at higher values of /3 does not necessarily 
make planet formation by GI possible in the inner regions of 
protoplanetary discs, since it is likely that /3 increases very 
rapidly towards the central star. The cooling time scale is 
proportional to (see Kratter et al.|20T0 l 

^ ^'''^ (95) 



^cool OC 
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where Kt is the opacity per unit mass. For an opacity law 
Kt oc et'', we get, using the ideal gas law: 



tc 



OC S 



(96) 



In a gravito-turbulent state, Q is constant, so we must have 
that E cx CaQ, and therefore: 



J. a + 26 — 4y-.3-t-a 



(97) 



In the outer re gions of protoplan etary discs, we expect a = 



and 6 = 2 (see [Bell fc Lin||l994[ ), which leads to /3 oc f2 



R, 



-9/2 



Since 13 is such a steep function of radius, it is very 



difficult to fragment at short distances, even in view of the 
results of this paper. However, because of possible inward 
migration, there is no need to fragment closer in, provided 
that the fragments can stay in the planetary mass regime. 
In principle, the mass of a fragment can even decrease if 



it migrates inward because of tidal stripping ( [Boley et al. 
2010} |Nayakshin|2010| ). 



Another way of saying that there is no critical cooling 
time scale for fragmentation in the simple case of /3-cooling 
studied here, is that there is no well-defined maximum stress 
the disc can sustain ( Rice et al.||2005 |. Formally, Qmax ~ 0. 
In practice, however, fragmentation becomes very rare for 
high values of /3, but in general great care has to be taken in 
using either /3c or Omax. There is no sharp boundary between 
discs that show fragmentation and discs that do not. 



8 CONCLUSIONS 

In this paper, we have studied the numerical convergence of 
the determination of the critical cooling time scale for disc 
fragmentation. We have seen that, in the two-dimensional 
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local approximation with a simple cooling law, there is 
no sharp boundary in terms of the cooling time scale be- 
tween discs that fragment and discs that do not. A 'steady', 
gravito-turbulent state consists of weak shocks as well as 
transient clumps, that will contract on a cooling time scale. 
If such a transient clump can survive in the turbulent back- 
ground for long enough, it will collapse and form a bound 
fragment. Since the weak shocks affect the disc only very 
locally, it is possible in principle for clumps to survive for 
many dynamical time scales. This makes disc fragmentation 
a stochastic process: most of these transient structures will 
be destroyed, but in the end one lucky clump will make it 
into a fragment. Transient clumps were found for all cool- 
ing times considered, and therefore fragmentation is possible 
in principle for cooling times up to /3 = 50. However, frag- 
mentation becomes increasingly rare for longer cooling time 
scales. 

It is possible that the efficiency of stochastic disc frag- 
mentation is affected by the approximations used. It remains 
to be seen whether similar effects can be observed in three- 
dimensional simulations with a more realistic cooling pre- 
scription. 
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